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Abstract: Seven next-to- leading order QCD evolution programs are compared. 
The deviations of the results due to different theoretical prescriptions for truncating 
the perturbative series are clarified, and a numerical agreement between five codes of 
better than 0.1% is achieved. Reference results for further comparison are provided. 

1 Introduction 

In order to exploit the full potential of HERA for deep-inelastic scattering (DIS), the highest 
possible luminosities and considerable efforts for the reduction of experimental systematic un- 
certainties are necessary. This will finally allow a measurement of the proton structure function 
F 2 over a wide range, with errors on the level of very few percent |TJ. To make full use of such 
results, and to allow even for combined analyses using the high-precision fixed-target data 
as well, the structure function evolution programs required for the necessary multi-parameter 
QCD fits have to be numerically and conceptually under control to a much higher accuracy. At 
least one order of magnitude is desirable. This accuracy is necessary to safely rule out contribu- 
tions to the theory error of a s (M|) which arise from the particular technical implementation of 
the solution of the NLO evolution equations. Due to the current apparent difference in a s (M§) 
as determined in e + e~ and DIS experiments ||, this question is of particular importance for 
the future QCD analyses based on the HERA structure function data. 

So far no high-precision comparison of next-to-leading-order (NLO) programs has been per- 
formed including the full HERA range. In previous studies partial comparisons were carried out 
demanding a considerably lower accuracy (see e.g. ref. ||). Other comparisons focussed on the 
valence range and compared the effect of different codes used for the QCD fit on Aq CD only, see 
ref. J4j. The required accuracy cannot be easily reached by just comparing results to published 
parametrizations, due to their inaccuracies, caused by respective numerical representations. 
Often also the physical and technical assumptions made are not fully documented. 

In this paper, we present the results of a dedicated effort, comparing the results of seven 
NLO codes under perfectly controlled conditions. The paper is organized as follows. In Sec- 
tion 2 we recall the basic formulae, and sketch the most commonly used approaches to the 
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evolution equations. Section 3 compares the differences of six 'global' evolution programs. The 
clarification of the deviations found there, using also a seventh program especially suited for 
the 'local' evolution of _F 2 , is described in Section 4. The size of the numerical differences 
which persist after this development is investigated in Section 5, where also reference results 
for further comparison are provided. Finally Section 6 contains our summary. 



2 Approaches to the Next-to-Leading Order Evolution 

The evolution equations for the parton distributions f(x, Q 2 ) of the proton are given by 
<9/(x,Q 2 ) 



<91nQ 2 



a s (Q 2 )P (x)+a 2 s (Q 2 )P l (x) + O(a s s ) <g> f(x,Q 2 ) . (1) 



Here x stands for the fractional momentum carried by the partons, and <g) denotes the Mellin 
convolution. For brevity, we have introduced a s (Q 2 ) = a s (Q 2 )/A7c. Eq. is understood to 
represent, in a generic manner, the non-singlet cases as well as the coupled quark and gluon 
evolutions. P^ and P^ 1 ' denote the corresponding leading order (LO) and NLO splitting 
functions, respectively (see, e.g., ref. ||). Only these two coefficients of the perturbative series 
are completely known so far, hence the solution of the evolution equations is presently possible 
only up to NLO. To this accuracy, the scale dependence of the strong coupling a s (Q 2 ) reads 

= -Poa 2 s (Q 2 ) - /^(Q 2 ) + 0(aj) . (2) 

Throughout our comparisons, we will identify the renormalization and factorization scales with 
Q 2 , as already indicated in eqs. (|l|) and (0). For different choices see refs. || |7j. Introducing 
the QCD scale parameter A, the solution of eq. (0) can be written as 

W ) -(3oHQ 2 /A 2 ) (31 In 2 (Q 2 / A 2 ) ' W 

Two approaches have been widely used for dealing with the integro-differential equations 
([]]). In many analyses, they have been numerically solved directly in x-space. We will exemplify 
some techniques applicable in this case for one particular program, choosing 'Qcdnum', which 
is based on the programs of ref. ||, and is planned to become publicly available |§. See, for 
example, ref. IIJ for a description of a differing a;-space implementation. 

In Qcdnum, the Q 2 evolution of the parton momentum densities is calculated on a grid in 
x and Q 2 , starting from the ^-dependence of these densities at a fixed reference scale Qq. The 
logarithmic slopes in Q 2 are calculated from eq. ([!]). To compute the convolution integrals, the 
assumption is made that the parton distributions can be linearly interpolated (at all Q 2 ) from 
one x gridpoint to the next. With this assumption the integrals can be evaluated as weighted 
sums. The weights, which are essentially integrals over the splitting functions, are numerically 
calculated (by Gauss integration) to high precision at program initialization. From the value 
of a given parton distribution and the slopes at Ql, the distribution can be calculated at the 
next gridpoint Q\ > Qq (or Q\ < Ql). This distribution then serves to calculate the slopes 
at Q\ etc., and the evolution is continued over the whole x-Q 2 grid. The evolution algorithm 
makes use of quadratic interpolation in \aQ 2 . 
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In this way, a fast evolution of parton densities is obtained, entirely based on look-up 
weight tables which are calculated at program initialization. The numerical accuracy depends 
on the density of the x grid and, to a lesser extent, on that of the Q 2 grid. In the comparisons 
presented here, 370 gridpoints in x covering 10 -5 < x < 1 have been used: 230 points distributed 
logarithmically for x < 0.2, and 140 points distributed linearly for x > 0.2. A logarithmic Q 2 
grid with 60 points covered the range 4 < Q 2 < 10 4 GeV 2 . 

An important alternative to the direct x-space treatment, employed in the analyses of refs. 
Ill], m based upon ref. fll3l , is to transform the evolution equations to Mellin-iV moments. 
The main virtue of this transformation is that the convolution is reduced to a simple product. 
Hence eq. ([!]) turns into a system of ordinary differential equations at fixed N, which allows 
for an analytic solution. Rewriting the evolution equations in terms of a s = a s (Q 2 ) using eq. 
(0), and expanding the resulting r.h.s. into a power series in a s , one arrives in NLO at 



d f( x > a ^ = L[p(o) f ^ + a (po-)( x )-0ip<p)(x 



0{a 2 s 



<8> f(x,a a 



(4) 



After transformation to iV-moments, its solution can be written down in a closed form for the 
non-singlet cases, with a = a s (Ql), as 



1 - 



is-aof p (i) 



r>y-> 
N 



Pi 

/V 



a£\ 
a J 



/jv(ao) • 



For the notationally more cumbersome, corresponding relation for the singlet evolution, 
reader is referred to refs. 



(5) 



the 
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From these analytic solutions, one can acquire the x-space results by one contour integral in 
the complex iV-plane, see ref. |l2j . Using a chain of Gauss quadratures, a numerical accuracy of 
this integration at better than 10~ 5 is readily achieved. In our comparisons, at most 136 fixed 
support points at complex iV-values have been used, with this maximal number employed only 
for very large values of x [R]. Due to the required non-trivial analytic continuations of the 
NLO anomalous dimensions [12|, this approach is technically somewhat more involved than the 
numerical x-sp&ce solution. On the other hand, since the Q 2 integration is done in one step, 
regardless of the evolution distance, and the use of fixed support points allows for performing 
the calculation of the anomalous dimensions only once at program initialization, this method 
is competitive in speed to the x-space iterations. 

A partly independent iV-space program has been developed during this workshop ||15|| , im- 
plementing an iterative numerical solution of the Mellin-transformed eq. (f|). Since one of the 
advantages of the iV-space approach is not exploited here, this program is so far not compet- 
itive in speed with the ones discussed before. It has however been of considerable value for 
cross-checks and theoretical investigations, see below. 

Before we now turn to the comparisons, it should be emphasized that a perfect agreement 
between the results based upon eqs. (p, (|), and (|5]) is not to be expected, since they all differ 
in terms of next-to-next-to leading order (NNLO), hidden under the O(a^) and 0(a 2 ) signs. 



3 The Initial Comparisons 

All our comparisons are performed under somewhat simplified, but sufficiently realistic condi- 
tions. We assume four massless flavours, in eq. (P as well as in eq. (0), at all scales considered, 
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i.e. effects due the non-zero charm mass and the existence of the bottom quark are not taken 
into account. All our results below will refer to the MS renormalization and factorization 
schemes, and the corresponding scales are identified with Q 2 . The reference scale Qq for the 
evolution, and the four-flavour QCD scale parameter A in eq. (D) are chosen as 



Ql = AGeV\ A^ = 250 MeV. (6) 
The following initial conditions are selected for the (anti-) quark and gluon densities: 



xu v (x, Qq) 
xS (x, Qq) 
xg(x,Q 2 ) 



A u x 
[xE 



0.5/ 



xUiyj x du j (^x j q J 



-0.2/ 



A d x°-°(1 



x 



= A s x-^{1 
A g x-°- 2 (l -x) 5 , xc(x,Ql) = xc(x,Ql) 



x 



(7) 



. 



The SU(3)-symmetric sea S is assumed to carry 15% of the nucleon's momentum at the input 
scale, and the remaining coefficients A{ are fixed by the usual sum rules. Finally Fi is determined 
by simply convoluting the resulting parton densities with the appropriate coefficient functions. 



Q 2 = 100 GeV 2 
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Figure 1: The differences between the up-valence, singlet quark and gluon densities, u v , E and g, 
and the proton structure functions F 2 , as obtained from evolving the input (0) with various NLO 
evolution programs |Tp|, |TJ, [15], [T5|, [T7[] to Q 2 = 100 GeV 2 . All results have been normalized to 
those of ref. 101. For a detailed discussion see the text. 
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The results of our first comparisons are shown in Figure 1. One notices the very good 
agreement between the programs @, used by the HERA collaborat ions. The differences 
are always much less than 1%, and the curves can hardly be distinguished, except for large 
x, with the present resolution. A similarly excellent agreement is seen between the two N- 
T5|j , except for very low x, where offsets up to 1.5% show up. The most 



space programs 

striking feature of the figure, however, is the very sizeable differences between these two groups 
of programs: the scaling violations, increasing (decreasing) the distributions at small (large) 
values of x, are considerably stronger in the results of refs. [14], 15| , although, of course, the 
same values for ot s are employed as in refs. || [TTJ. This effect reaches a magnitude of as much 
as 8% for the structure function F 2 at the smallest x-values considered. 

As stated above, perfect agreement had not been expected due to theoretical differences, 
but the size of this offset was a surprise to most of us. It initiated quite some checking and 
programming activity, which will be summarized in the next section. 

Also shown in the figure are the results obtained by the x-space evolution programs of the 
MRS and CTEQ global fit collaborations fy| |17| . Very good agreement to the results of refs. 
|| 10] is found for the valence quarks, except for ref. [16| at extremely low values of x. In the 



singlet sector, however, significant differences are observed for some quantities: 1.5 - 3% on the 



gluon density in ref. |16l, and up to 4% on the sea quark distributions in ref. |L7| 



4 Pinning Down the Differences 



Besides checks and comparisons of the numerical values of the NLO splitting between the codes 
of refs. || [T3|, [L5| , the large differences discussed in the previous section led to three program 
developments, which together allowed for their full understanding on an unprecedented level. 



A program for a local representation of the evolution of F 2 close to the initial scale, 
completely independent of all previous ones, was added to the comparisons 



The code of ref. (jOJ was extended to include, still in moment space, an option for evolving 
also on the basis of eq. (|TJ) instead of (|). 



The program of ref. JT4J was used to simulate an iterative solution of eq. (f|) as performed 
in ref. ]T5| , and additionally two new iterative options, one of them equivalent to eq. (^), 
were introduced into this package. 



The results of these efforts are displayed in Figure 3, where we show the evolution of F 2 , close 
to our reference scale Ql = 4 GeV 2 , for three typical values of x. The differences, depicted 
in the previous figure for Q 2 = 100 GeV 2 , build up very quickly near Q 2 ,: already around 
10 GeV 2 they are close to their final level. The second important observation is the perfect 

which immediately 
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agreement of the local representation [18] with the x-space codes [E 
stopped any speculations on possible problems in the latter programs. Next one notices that 
the 1% small-x difference between refs. []14j and [|15j is perfectly understood in terms of the 
slightly different contributions truncated away in eqs. (f|) and (||), cf. ref. ||. The concluding 
step is the comparison of the modified evolution of |T5[ with the results of || fTDj, [TJ|. This 
reveals that in fact virtually all offsets between the results of refs. [EJ, 1C, 14, 15] in Figure 1 are 
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Figure 2: A comparison of results on the Q 2 evolution of F 2 close to the reference scale Ql = 4 GeV 2 . 
The results of refs. [01 O, |I2| (denoted by B, R, and V) are as in the previous figure. BvN 
represents a local representation of the F 2 evolution [|18[], and the curves R' and V check the 
numerical consistency by adapting to the theoretical assumptions of refs. and |15|| , respectively. 



due to the differences introduces by the employed truncation prescriptions for the perturbative 
series the NLO level, i.e. by terms of NNLO and beyond. 



The origin of the differences between the results of refs. [16|, [L7| and our programs could 
not be clarified during this workshop. Hence for the very precise comparisons to which we now 
turn, we will keep only our five program packages, which agree, at least, sizeably better than 
to 1%. 



5 The Achieved Numerical Accuracy 

Armed now with at least two different codes for any of the truncation prescriptions of Section 2, 
we can proceed to explore the limits of the agreement of our five program packages under 
consideration. This complete coverage will be used for comparing all programs, even those 
with conflicting theoretical treatments, in one figure. For this purpose, the results of refs. 
T0| have been normalized to the modified evolution of ref. []T3] (based upon eq. (0)), whereas 



the 'iterated' evolution of ref. |14| is normalized to the original results of ref. |15] (based upon 

The results are shown at two fixed Q 2 values in Figure 3 for the parton distributions, and 
in Figure 4 for Q 2 evolution of the proton structure function F 2 at three fixed values of x. The 
total spread of the results at Q 2 = 100 GeV 2 amounts to at most about 0.05%, except for very 
large x, where the distributions, especially the gluon density, become very small. Even after 
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evolution to 10 4 GeV 2 , the differences are still on the level of 0.1%, meeting the goal formulated 
in the introduction. Moreover, there is no reason for failing to reach an even higher accuracy, 
at least to 0.02% as already achieved between the iV-space programs, also in x-space, e.g. by 
increasing the still not too high number of Q 2 grid points in the program of ref. |J. 

Finally, for the convenience of those readers who want to check their own existing of forth- 
coming NLO evolution program to an accuracy well below 0.1% over a wide range in x, we show 
in Table 1 two sets of reference results, which represent the evolution of the initial distributions 
(0) under the conditions (^), according to eq. (0) and eq. (|5|) to Q 2 = 100 GeV 2 . 
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Table 1. Reference results at Q 2 = 100 GeV 2 for the NLO evolution using the direct solution of 
eq. (1) (upper half), and the truncated analytic solution (5) (lower half). The initial conditions are 
specified in eqs. (6) and (7). The estimated numerical accuracy of these results is about 0.02%. 



6 Summary 

The results of seven programs for the NLO evolution of parton densities and structure functions 
have been compared. Differences due to terms of NNLO, truncated differently in the various 
implementations, turn out to be larger than anticipated. They can reach, e.g., about 6% at 
x = 10~ 4 and Q 2 = 100 GeV 2 . A full quantitative understanding of these differences has been 
achieved at an unprecedented level of accuracy for five of these codes. There the remaining 
numerical differences are on the level of ±0.02% at Q 2 = 100 GeV 2 . Two sets of reference 
results, according to different theoretical prescriptions, have been provided for further high- 
precision checks of evolution programs. 
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Figure 3: The remaining numerical ^-dependent deviations on u v , £ and g at Q 2 = 100 and 10 4 
GeV 2 between the programs of refs. [[§, 0, [14], |15j, after the differing theoretical assumptions have 
been corrected for, see the text. The results have been normalized to those of ref. JT5| . 
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Figure 4: The residual relative offsets in the Q 2 evolution of F 2 between the programs of refs. 



14] , |15fl , after removing the effects due to the different theoretical NLO prescriptions, for three 



typical values of x. As in all previous figures, the initial distributions are taken from eq. (0). 
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